Solar Physics 
DOI: 10.1007/' 



A Bayesian Analysis of the Correlations Among 
Sunspot Cycles 

Y. Yu 1 • D.A. van Dyk 2 • V.L. Kashyap 3 • 
C.A. Young 4 

s ! ' © Springer •••• 

o ■ 

^v^j . Abstract Sunspot numbers form a comprehensive, long-duration proxy of so- 

lar activity and have been used numerous times to empirically investigate the 
properties of the solar cycle. A number of correlations have been discovered over 
the 24 cycles for which observational records are available. Here we carry out 
a sophisticated statistical analysis of the sunspot record that reaffirms these 
correlations, and sets up an empirical predictive framework for future cycles. An 
advantage of our approach is that it allows for rigorous assessment of both the 
fV^ statistical significance of various cycle features and the uncertainty associated 

. with predictions. We summarize the data into three sequential relations that 

estimate the amplitude, duration, and time of rise to maximum for any cycle, 
O I given the values from the previous cycle. We find that there is no indication 

I ' of a persistence in predictive power beyond one cycle, and conclude that the 

^ . dynamo does not retain memory beyond one cycle. Based on sunspot records 

" up to October 2011, we obtain, for Cycle 24, an estimated maximum smoothed 

Ci • monthly sunspot number of 97 ± 15, to occur in January-February 2014 ± 6 

months. 

>: 

^sO . 1. Introduction 

°: 

, . Sunspot numbers constitute the longest continuous record of observations in 

astronomy, having been recorded in observatories worldwide since the Dalton 
minimum, and they are available as monthly estimates since 1749. Because of the 
' multi-generational span, both the quality of the observations and the techniques 

— , . used to record them have varied, and thus these data present a challenge for 

interpretation. The data are maintained at the Solar Influences Data Analysis 
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Center in Belgium (http://sidc.oma.be) as the international sunspot numbers 



(SSNs) and their suitability as a proxy for solar activity has been found to be 
good by comparison with other proxies (Waldmeier, 1971). (Recent efforts to 
recalibrate the SSN have suggested potential differences between the different 
proxies (e.g. Svalgaard, 2010). However, we do not include these corrections in 
our analysis because first, there is a possibility that a fundamental change has 
occurred in the manifestation of the activity proxies in recent years, and second, 
there is no evidence that the calibration needs to be changed over the entire 
historical record.) The long duration of the dataset allows us to analyze a much 
larger number of cycles than with other, more recently developed proxies of solar 
activity (such as sunspot area, umbral fields, 10.7 cm flux, etc.). 

Accurate prediction of solar activity cycles is an important area of research 
since variations in "space weather" caused by solar activity affect radio commu- 
nication, the performance of low-Earth orbit satellites, and geomagnetic activity 
(e.g. aurorae). Studying the behavior of sunspot cycles is important not only for 
understanding the physics of solar activity, but also for the planning of space 
missions. Sunspot numbers are valuable as an indicator of solar activity, and 
identifying recurring patterns in sunspot cycles is therefore crucial from both an 
empirical and a theoretical perspective. Indeed, the 11-year activity cycle of the 
Sun was first discovered by noticing the same cycle in SSNs (Wolf, 1852). While 
correlations with solar activity have been identified in other indicators (solar 
flare numbers, sunspot areas, 10.7cm flux, etc.; see Hudson, 2007) and studies of 
solar activity have been extended back many millenia using dendrochronological 
data (Bonev et al, 2003; Solanki et al, 2004), the SSN data are the first rung 
in the ladder to calibrate all observations of the solar activity cycle. Waldmeier 
(1935, 1939) noted that sunspot cycles tend to take less time to rise to maximum 
than to fall to minimum. Other important relations, such as the correlation 
between the duration of a cycle and the amplitude of the next cycle (amplitude- 
period effect), can be used to predict characteristics of the upcoming cycle years 
in advance (e.g. Hathaway et al., 1994, 2002; Watari, 2009). Here we analyze 
the SSN data to derive statistically meaningful phenomenological correlations 
between the various parameters that observationally define a solar cycle. Such 
correlations act as constraints on theoretical models of dynamo action that seek 
to explain solar activity (Schiissler, 2007). Analysis of long-duration activity 
cycles (Usoskin et at, 2007) suggest that they are driven by a stochastic or 
chaotic process, and it is therefore necessary to identify the statistical properties 
of the activity cycles in the current era. 

Prediction methods for upcoming cycles include those based on i) solar dy- 
namo models (Choudhuri, 1992; Charbonneau and Dikpati, 2000; Dikpati et 
al., 2006; Choudhuri et al., 2007; Charbonneau, 2007), ii) precursors such as 
geomagnetic aa indices (Hathaway and Wilson, 2006), and iii) statistical analyses 
with historical data (Hathaway et al, 1994; Benestad, 2005; Xu et al, 2008; Gil- 
Alana, 2009), among others. Pesnell (2008) gives a review of a large range of 
predictions made for the upcoming Cycle 24. Given the current debate over the 
amplitude of Cycle 24, for which different physical models yield substantially dif- 
ferent predictions (see, for example, Dikpati and Gilman, 2006, and Choudhuri et 
al, 2007), a statistical method that uses only the SSN data will provide a useful 
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benchmark for comparison. Methods based on statistical extrapolation, however, 
rely on various assumptions just as physical models do, and may not fully account 
for the uncertainties involved. For example, for Cycle 23, the smoothed maximum 
sunspot number as predicted by several researchers varied considerably, from 80 
to 210. Kane (2001) noted that among 20 predictions, only 8 were within a 
reasonable range of the actual value. 

We describe our analysis of SSN data in Section [51 and in particular describe 
the statistical model used to determine the correlations between the parameters 
defining a given cycle from the previous cycle (Section 12. ip and fitting the model 
to the data (Section 12. 2[) . We report the resulting correlations and discuss their 
implications in Section [3J and summarize our work in Section 2) 

2. Analysis 

2.1. Two-Stage Statistical Model 

This article proposes a two-stage statistical model that accounts for the uncer- 
tainty both in smoothed monthly sunspot numbers and in predicting future cycle 
characteristics from historical data. In the first stage of the statistical model, 
cycle characteristics, such as amplitude, duration, and rising time, are estimated 
from raw SSNs. Then, in the second stage, relations between the characteristics 
of consecutive cycles are examined. This results in three sequential relations 
that summarize known features of sunspot cycles such as the Waldmeier effect 
and the correlation between amplitudes of successive cycles. These statistical 
properties place a constraint on any physical model that attempts to explain the 
solar cycle behavior. 

In this section we fit the two stages of the statistical model separately. This 
involves first modeling the cycles in Section 12.1.11 and then modeling the rela- 
tionships among the parameters that describe the cycles in Section 12.1.21 We 
discuss how these two stages can be combined into a single coherent statistical 
model in Section 12.21 and show the results of this coherent fit in Section [31 

2.1.1. Stage One: Modeling the Cycles 

When extracting cycle length, rising time, and amplitude information from SSN 
data, we adopt an approach similar to the two-parameter curve fitting of Hath- 
away et al. (1994; see also Sabarinath and Anilkumar, 2008, and Volobuev, 2009, 
who propose other functional forms). For cycle i, suppose is the starting time, 
£max is the time of the cycle maximum, is the end time, Cj is the amplitude, 
and Ut is a parameter that captures the "average solar activity level" at time t. 
We postulate that 

(i) 

• for the rising phase t < t lnax 




and 



(1) 
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for the declining phase t > t 



(0 



u t = a[i- 



t - 1 



(0 



(2) 



where a±,a2 > 1 are shape parameters assumed to be the same for all cycles. 
A curve described by the two postulates is illustrated in Figure [TJ We do not 
assume that the starting point of the next cycle is identical to the end 

point of the current cycle, t\ . When tq < t\ , the two cycles overlap, and 
the activity level [Ut] during the overlapping period is defined as the sum of 
the contributions of the form ([2]) and ([I]) from these two cycles. We adopt this 
parameterization because it is simple, flexible, easy to interpret, and fits the data 
well. 




(i) 



2040 



2042 



2044 2046 
year 



2048 



2050 



Figure 1. Parameterized form of a solar cycle. We illustrate (7 t 2 with c; = 10, a± = 1.9 
and ct2 = 1-1) where Ut is specified by {TJ and {2}- Because we model after a square-root 
transformation, Uf approximates the shape of a cycle on the original scale of the SSN data. 

Given a total of k = 25 cycles (including the incomplete Cycles and 24), 
we relate the parameters a — {a\,a2) as well as To = {to \ i = 0, — 1), 
T max = (tma X , i = 0, . . . , k - 1), Ti = (ti' , i = 0, . . . , k — 1), and C = (ci, i = 
0, . . . , k — 1), to the observed data by a linear model: 

y/Yi\(To,T mm ,T 1 ,C,a,p i o 2 ) , & N(fi + U t , a 2 ), (3) 

where the parameter /3 may be regarded as a baseline, and we model Y t , the 
monthly average sunspot number at time t, after a square-root transforma- 
tion. This transformation is used to stabilize the variance, since higher sunspot 
numbers are also associated with higher variability. (A variance-stabilizing trans- 
formation is a mapping [/(y)] of the data [y] such that the variability of f(y) 
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is constant relative to its mean value. For example, for counts data that follow 
a Poisson(/i) distribution, the variance increases with the mean /i. The variance 
of y/y, however, is approximately the same for different /i. The situation with 
SSNs is similar and a square-root transformation is suitable. Note that results 
such as predictions are easily obtained on the original scale by inverting the 
transformation. ) 

Our approach differs from other curve-fitting methods (e.g. Hathaway et al., 
f994) in that we model all cycles jointly and we estimate the starting and 
ending points of the cycles from the data rather than fixing them in advance. In 
addition, we model the monthly sunspot numbers directly rather than smooth 
them beforehand, e.g. by a moving average, as is often done in other statistically 
based prediction methods. It may be interesting to analyze daily data in a future 
work; see Noble and Wheatland (2012) for an analysis of the daily fluctuations. 
One intrinsic difficulty with daily data, however, is that the same spot or spot 
group is counted every day until it disappears or rotates away over the limb 
(aside from the issue of merging or splitting spots). We choose the monthly data 
partly to alleviate this problem. 

When To and T max are treated as unknown parameters, model ([3]) is not the 
usual linear regression model that can be easily fit by ordinary least squares. To 
overcome such difficulties, we adopt a Bayesian approach and use Markov chain 
Monte Carlo (MCMC) to simulate samples from the posterior distribution. We 
discuss Bayesian analysis, MCMC, and how we use them to fit our model in 
Section 12.21 The fitted model and residuals are illustrated in Figure [5J These 
are based on a random draw from the posterior distribution after fitting ©. 
Plots using multiple posterior draws are qualitatively the same and are omitted. 
The residuals are defined as \/Y~t — j3 — Ut as in ((3J) . The residual plot reveals a 
reasonably good fit, given the simple functional form of (JTJ) and fl2J). 

Note that (J3j> alone docs not specify relationships between consecutive cycles, 
and is therefore not useful for predicting cycle characteristics of entirely new 
cycles. In Section 12.1.21 we explore relationships between cycle characteristics, 
and build additional structure on Equation (|3|) to enhance its usefulness in 
predictions. 

Table [1] displays the fitted values and standard errors for some parameters 
that are common to all cycles. The difference between a,\ and a 2 suggests an 
asymmetry in shape between the rising and declining phases of a cycle. We have 
also fitted with individual (ai, a^) for each cycle, but the results are similar and 
hence omitted. 



Table 1. Fitted values (posterior means) and 
standard errors (posterior standard deviations) for 
ai, «2, and cr 2 . 



a l 


"2 


P 


1.9 ±0.11 


1.1 ± 0.04 


0.47 ±0.18 1.2 ±0.03 
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SSNs with fitted values 
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year 
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Figure 2. Fitted values and residuals for the regression model in Q for one posterior sample. 
Top: sunspot numbers with fitted values of T max represented as vertical lines (only the last few 
cycles are shown for illustration). Bottom: histogram of residuals with approximating normal 
density curve. 

Table [2] summarizes the parametric fits to each cycle. (Again, the model-fitting 
procedure will be described in Section l2~2"n Note that the interval between when 
the old cycle ends and the new cycle begins is usually negative, suggesting that 
the new cycle begins before the old cycle ends. This is consistent with recent 
results based on torsional oscillations (Hill et al., 20f0). 

2.1.2. Stage Two: Modeling Relationships Between Consecutive Cycles 

The goal at this stage of the statistical model is to build up an empirical mech- 
anism that generates the amplitude, duration, and rising time of a cycle, given 
those of the previous cycle. As a preliminary check, we carry out correlations 
between different parameters that define the model in a given cycle as well as 
between adjacent cycles. The results of the correlation analyses arc reported 
in Tables [3][5l The tables list Spearman's rank coefficient [p] (Kendall, 1975) 
computed using the best-fit parameter values and the corresponding p-value. 
(The error bars on p are computed from 200 posterior samples of the parame- 
ters and represent the robustness of the correlation. The p- value represents the 
probability that such a correlation may be obtained by chance.) Within a cycle 
(Table a clear anti-correlation is seen between the amplitude [cj] and the 
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Table 2. Cycle-profile parameter posterior means for each cycle 



cycle amplitude rise [yr] fall [yr] gap [yr] start [yr] end [yr] 





Ci ' 


- VSSN 


tmax ^0 


tl ^inax 


t 


d ~ h 




to 




ti 





9 


_ _ _i_n 


36 
34 


5 


00±|;|1 


6.43+81° 

— 0.35 


-1 


■40+811 

— 0.18 


1745 


0+2.6 

— 3.2 


1756. 


fi +0.2 
— 0.2 


1 


8, 


- -4-0 

.55+°; 


2:5 
22 


6 


. „ _i_0 1 9 
/ti -0.32 


« 2 o + 0.44 
" —0.40 


-2 


in+0.26 
— 0.24 


1755 


9+0.2 
— 0.2 


1768. 


q+0.3 
' —0.3 


2 


10 


321°, 


23 
21 


3 


^ x -0.25 


7 no+0.37 

/ . OO n ni 

—0.31 


-1 


■io+815 

— 0.24 


1766 


o+O.l 
—0.1 


1777. 


5 +0.3 

— 0.3 


3 


12 


,-..--1-0 

.251°; 


2r, 
24 


2 


„ ._i_0 21 

u ^-0.2i 


7 4 o+0. 26 
" —0.24 


-1 


,07+°'}6 

—0.18 


1776 


.4+81 

—0.1 


1786. 


o+°-? 

— 0.2 


4 


11 


.571°' 


24 
22 


3 


U '-0.24 


11.73+811 

— 0.24 


-1 


1 c + 0.24 
—0.18 


1784 


9+8'? 

— 0.2 


1799. 


7 +0.2 

— 0.2 


5 


6 




2') 
24 


5 


yy -0.32 


5.59_q ^2 


1 


CT2 + - 48 

— 0.44 


1798 


6+81 
—0.1 


1810. 


i+0.2 
— 0.2 


6 


6 


,20+° 


2~> 
22 


5 


no i vj ."±u 
yo -0.45 


5 . 88"^"n •> 

—0.38 





Q7+0.26 
—0.32 


1811 


7 +0.3 
—0.2 


1823. 


cr + 0.3 
— 0.2 


7 


8, 


^ « -t- n 
33+°; 


24 
22 


6 


° z -0.23 


4.73+81! 

— 0.24 


-0 


7 o+0.14 
1 —0.19 


1823 


c+0.1 
,u — 0.1 


1835. 


1+0.2 
— 0.2 


8 


11 


631°, 


27 
2- r > 


2 


81 -o!l6 


8.43+81? 

—0.35 


-1 


—0.16 


1834 


.4+81 

—0.1 


1845. 


6+81 
— 0.3 


9 


10 


a n 4-0 

48l°' 


23 
24 


1 


„~4-Q 25 
yz -0.25 


71 +0.37 
—0.30 


-1 


,00+°'?| 

—0.17 


1843 


9+0.2 
— 0.2 


1857. 


g+0.3 

— 0.2 


10 


9 


- -1-0 

,811°; 


2~i 
21 


3 


„„_|_0 94 


9.16+81? 

— 0.41 


-1 


.79+8'?° 

— 0.21 


1856 


,6+81 

— 0.2 


1869. 


£+0.3 

— 0.2 


11 


11 


05+° 

— 


2.') 

2(j 


3 


2n +0.13 
— 0.20 


8.07+8 -?S 

— 0.24 


-0 


9O + 0.19 
1 —0.22 


1867 


7 +0.1 
—0.1 


1879. 


0+8I 

— 0.1 


12 


8. 


231°' 


24 
22 


4 


79 +0.44 
' ' -0.47 


7.29l°;f 4 


-0 


60 +0.27 
,DU -0.23 


1878 


7+0.1 
'-0.1 


1890. 


7+0.3 
'-0.3 


13 


9 


151°, 


2:', 
2:s 


3 


47+0.28 
4 '-0.30 




-0 


fic; +0.23 
DO -0.19 


1890 


1+0.1 
!-0.1 


1902. 


6+°' 2 
°-0.2 


14 


8, 


ii+O. 
1L -0 


2(i 
25 


4 


Q o+0.27 
y °-0.40 


fi 77 +0.48 
D -' '-0.45 


-0 


ni+0.26 
ui -0.24 


1901 


O+0.2 
y -0.2 


1913 


7+0.2 
'-0.3 


15 


9 


.441° 


2(i 
2.'! 


4 


45+0-22 
^ d -0.28 


6 27+ ' 39 


-0 


05+0.12 

,aJ -0.21 


1913 


7+0.1 
' -0.2 


1924. 


4+0-3 
-0.2 


16 


8 


.821°' 


2!) 
2") 


4 


gi+0.39 
O1 -0.36 


6 7Q+ - 46 
o -' 3 -0.38 


-0 


77+0.19 
' ' ' -0.23 


1923 


4+0-2 
-0.2 


1934. 


o+0.2 
8 -0.2 


17 


10 


701° 


24 
2.'! 


4 


9C+0.25 
zo -0.17 


7 4 o+0.32 
'• 4c> -0.35 


-1 


97+O.I8 
z '-0.23 


1934 


n +o.i 

u -0.1 


1945. 


7+0.2 
'-0.2 


18 


12 


.491° 


27 

28 


4 


O9+0.23 
u -0.19 


6 73+ - 19 
°- '°-0.23 


-0 


51+0.17 
d± -0.08 


1944 


5+0.1 
°-0.1 


1955. 


9 + 0.1 
-0.1 


19 


14 


.191° 


2!) 
2(1 


3 


QO+0.10 

zo -0.15 


7 81+ ' 19 
'• o± -0.16 


-1 


Q4+0.12 
U4 -0.13 


1954 


7+0.1 
'-0.1 


1965. 


0+0.2 

8 -0.2 


20 


10 




21. 
24 


1 


oq+0.29 
^"-O.SO 


Q Q7+0-53 

y-y '-0.47 


-2 


nfi +0.31 
UD -0.29 


1964 


7+0.1 
'-0.1 


1979. 


1+0.3 

-■■-0.4 


21 


13 


.041° 


2 r > 
2<) 


3 


46 +0.21 
^°-0.21 


7 09 + O.26 
-0.32 


-1 


m+0.09 

ul -0.16 


1977 


+o.i 

u -0.1 


1987. 


o+0.2 
S -0.2 


22 


12 


.911° 


2!) 
24 


3 


49+0.18 

-0.25 


7 1O+0.22 
' - la -0.28 


-0 


00+0.13 
s»_0.12 


1986 


0+0.1 

8 -0.1 


1997. 


4+0.2 

^-0.1 


23 


10 


.901° 


2(i 
26 


4 


4u -0.23 


01 +0.36 
°- OL -0.39 


-0 


90+0.34 

ZD -0.33 


1996 


c+0.2 
°-0.2 


2009. 


9+0.2 
z -0.2 



rise time [i^ax — t l ] , as well as between rise time and fall time [t\ — tj nax ] ; i. e. 
strong cycles rise to maximum quickly (Hathaway et ai, 1994), and when they 
rise quickly they tend to decline slowly. Correlations also appear to be present 
between amplitude and fall time, with strong cycles correlated with long declines 
(consistent with the above anticorrelations) , and inversely between fall time and 
cycle gap [<q +1 — 4]- However, the direct correlation between the rise time and 
cycle gap is not statistically significant. 

We show the correlations of a parameter with the parameter values of the 
following cycle in Table H] and with the parameter values in the preceding cycle 
in Tablc[5] Only the amplitude is significantly correlated across the cycles. There 
is weak evidence for a correlation between the current cycle amplitude and the 
fall time of the next cycle, and inversely between the previous cycle gap and the 
current amplitude. For more recent determinations of correlations across cycles, 
see Kane (2008), Kakad (2011), Ramesh and Lakshmi (2012); see also Vaquero 
and Trigo (2008) for a cautionary note. 

We thus consider the relation between parameters in consecutive cycles so 
that they may be exploited to predict sunspot numbers of an entirely new cycle. 
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Table 3. Correlations within cycle 





risetime 

fi _ fi 
''max ''O 


falltime 

fi _ fi 

L l ^max 


gap 


amplitude 

Ci 


p = -0.81 ± 0.06 
p = 0.00 


p = 0.48 ± 0.07 
p = 0.02 


p = -0.32 ± 0.06 
p = 0.13 


risetime 
t* — t i 

L max ''O 




p= -0.62 ±0.09 
p = 0.00 


p = 0.29 ± 0.09 
p = 0.18 


falltime 
fi _ fi 

L l ^max 






p = -0.46 ± 0.08 
p = 0.02 



Table 4. Correlations with following cycle 







amplitude+ 


risctimc+ 


falltimc+ 


gap± 


amplitude 


P 


0.46 ± 0.02 


-0.21 ± 0.06 


0.39 ± 0.05 


-0.06 ±0.07 




P 


0.03 


0.34 


0.07 


0.80 


risetime 


P 




0.05 ± 0.10 


-0.27 ± 0.07 


-0.06 ±0.13 




P 




0.83 


0.21 


0.37 


falltime 


P 






0.20 ±0.07 


0.19 ±0.07 




P 






0.36 


0.39 


gap 


P 








0.25 ±0.10 




P 








0.25 



Table 5. Correlations with previous cycle 







amplitude 


risetime 


falltime gap 


amplitude— 


P 


0.46 ±0.02 








V 


0.03 






risetime— 


P 


-0.17 ± 0.08 


0.05 ±0.10 






P 


0.44 


0.83 




falltime— 


P 


0.07 ± 0.05 


-0.08 ±0.06 


0.20 ±0.07 




P 


0.76 


0.72 


0.36 


gap- 


P 


-0.37 ± 0.07 


0.29 ±0.07 


-0.25 ± 0.09 0.25 ±0.10 




V 


0.08 


0.18 


0.25 0.25 



Wc have carried out numerous checks of the parameter combinations. In Figure [3] 
we display two such relations between the fitted values (posterior means) based 
on the first-stage model (Equation the amplitude of the next cycle [q+i] 
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against the amplitude of the current cycle [cj] (left plot), and against the time 
from maximum to the start of the next cycle [tQ +1 ^ — tmax] (right plot). 




Figure 3. Relationships between the next amplitude c;_|_i and the current amplitude c;, and 
between Ci_|_i and t — tm a x- Reported are posterior means based on the first-stage model 
in Equation (0. 

Based on the above exploratory analyses, we consider a combination of param- 
eters to define a statistical model for predicting the next cycle parameters. First, 
we enhance the predictive power of the positive correlation between successive 
amplitudes by combining it with the negative correlation shown on the right plot 
of Figure [3l 

c i+ i~ft+7i (i) +N(0,af). (4) 

Next, we observe the following form of Waldmeier effect relating the rising 
time of each cycle to its amplitude, 

t&x 1 * - ~ S 2 + 72 c i+1 + N(0, a 2 2 ). (5) 

This relationship is illustrated in the middle panel in Figure 01 again based on 
the first-stage model in Equation ([3]). 

Finally, we incorporate the relatively weak correlation between the amplitude 
and the duration of the declining phase of each cycle, 

- t^> ~ 6 3 + 73 c t+1 + N(0, at). (6) 

Figure H] illustrates Equations flU, © and ©. Although the third panel of 
Figure H] exhibits a relatively weak relationship, we note that Equation (jSJ) is 
not needed for predicting the timing and amplitude of the next solar maximum, 
which is often the main goal. 



SOLA: ssn_solphy_2012.tex; 9 August 2012; 0:17; p. 9 



Y. Yu et al. 



Note that only one of these relations explicitly connects the parameters of 
one cycle with that of the next (Equation (U)). The other two (Equations (O 
and (O) are based on parameter correlations within a cycle, but are computed 
for the following cycle. This sequence, of first computing the amplitude of the 
following cycle, and using that to compute the rise and fall times of that cycle, 
is an important facet of the predictive capacity of our model. The calculations 
cannot be carried out in a different order. 




Figure 4. Three linear relationships useful for prediction, based on fitting the first-stage 
model in Equation ((3J. 



It should be emphasized that, because of the large amount of raw SSN ob- 
servations, the quantities c*, imlx, tf, i = 0,...,23, and , i = 1,...,24, 
are well constrained by the first stage model in Equation ((3]). Nonetheless, their 
fitted values (posterior means) do not account for the uncertainties and arc used 
only for the purpose of illustration. The mathematical forms of the relationships 
in Equations Q - © are found by examining the posterior means of these 
parameters, but ultimately the parameters will be fit using the raw data. In 
Section [2.2l we describe a two-level joint modeling approach which automatically 
accounts for uncertainty in these quantities. 

We have also explored possible correlations between separated cycles (i.e., 
between the k th and the (/c±2) th cycles). However, given the characteristics of 
Cycle k — 1, we find no evidence for a conclusive dependence of Cycle k on Cycle 
fc±2. Such dependences are weak and have little predictive value. Therefore we 
focus on lag-one dependence. 
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2.2. Model Fitting 

We adopt a Bayesian approach based on the relations found in the previous stage 
and carry out principled parameter estimations for the correlation coefficients. 
We employ Markov chain Monte Carlo methods that fit both stages of the sta- 
tistical model simultaneously. This approach makes it easy to correctly account 
for the uncertainty both in estimating the parameters of the individual cycles 
and in modeling the relations among cycles. 

2.2.1. Bayesian Hierarchical Models 

The results in Section I2TT1 are obtained by fitting Equation ^ using a Bayesian 
approach. In this section we describe the model fitting procedure as well as 
how the Bayesian approach can be used to fit both stages (Equations © and 
(HI)"®) together as a single more coherent statistical model. See Gelman et al. 
(2004) for an introduction to Bayesian hierarchical modeling and the associated 
computational strategies. See van Dyk et al. (2001) for the use of Bayesian 
methods in the context of highly structured models for spectral analysis in 
high energy astrophysics, and Esch et al. (2004) in the context of multiscale 
image reconstruction. In a Bayesian analysis, the likelihood function [p(y|6>)] is 
combined with a prior distribution [p{9)\ to form a posterior distribution, 

p(6\Y)<xp(Y\e)p(6), 

and all inference is derived from this posterior distribution. One may build 
further structures on the prior by introducing some hyper-parameter [rj] with 
its own prior [p(i])] (the hyper-prior), and replacing p(9) with p(9\ri). This leads 
to a two-level model. By Bayes' Theorem, the joint distribution of (9, rj) given 
the data Y can be written as 

p(e, V \Y)xp(Y\9,ri)p(9\ri)p(ri), 

where p(Y\6, rj) is the likelihood of observed data (the first stage), p(9\rf) is the 
second stage distribution for 9, and p(rj) is the hyper-prior. Inference concern- 
ing 9, for example, is based on its marginal posterior distribution p(9\Y) = 
fp(e,ri\Y)dt]. 

In our case Y — {Y t } and 9 is the collection of parameters 

o = (ro,r max ,ri,c , ,a,/3,<7 2 ). 

The likelihood function p(Y\6) is determined by Equation ([3]). We impose in- 
dependent uniform prior distributions on j3 and logcr, as is commonly done for 
such parameters in a first-stage regression model. For each of a\, 0.2, a uniform 
prior on [1,3] is used to allow for a flexible range of cycle shapes. When fitting 
the first stage model ([3]) alone, as in Section [2. 1.11 noninformative uniform priors 
are assigned to the other components of 6, i.e. T ,T max ,Ti, and C, subject to 
natural constraints on their ranges. In this section we fit the two stages jointly, 
and these components of 9 are linked together by Equations ([!]) - © . In addition, 
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we express the starting point of the next cycle t^ +1 \ given the end point of the 
current cycle , as 

4*+!) ~t«+iV(0,r 2 ). (7) 

The parameter t 2 regulates how far apart ^ and t± are allowed to be. Allow- 
ing 4* to be different from offers additional flexibility. In summary, Equa- 
tions (|l])-([7} specify the distribution of Tb,T max ,Ti, and C given the hyper- 
parameters rj = (r 2 ,jj,Sj,<jj , j = 1,2,3). For these hyper-parameters we use 
independent non-informative priors (specifically, uniform priors) on T,~/j,6j,(Tj, 
j = 1,2,3. 

2.2.2. Model Fitting with MCMC 

Markov chain Monte Carlo (MCMC) techniques are used to draw samples from 
the posterior distribution, which can then be summarized as point estimates and 
error bars of the parameters of interest {e.g. van Dyk et al., 2001 and Park et al., 
2008). In general, if we can simulate random samples 0W, . . . , 9^ from p(9\Y), 
the target posterior distribution, then inferences for quantities of interest can be 
derived by examining the empirical distribution of these samples. For example, 
suppose 9 is one-dimensional and has a symmetric and unimodal posterior dis- 
tribution, then a natural point estimate of 9 is the posterior mean, which can be 
approximated by L _1 X)z^i ^ i th° empirical average of the posterior sample. 
The posterior standard deviation serves as a natural one-sigma error bar, and 
can be approximated by the sample standard deviation of 9^\ . . . , 9^ L '. 

In high-dimensional situations when direct simulation from the posterior dis- 
tribution is difficult, as is the case of our analysis, one may adopt an MCMC 
approach, which constructs a Markov chain with the desired posterior distri- 
bution p(9\Y) as its stationary distribution. After the Markov chain reaches 
equilibrium, the iterations of 9 can then be used as (dependent) samples from 
the target distribution. Two well-known methods for constructing such Markov 
chains are the Metropolis-Hastings (M-H) algorithm (Metropolis et al., 1953; 
Hastings, 1970) and the Gibbs sampler (Geman and Gcman, 1984; Gelfand and 
Smith, 1990). In an M-H strategy, to obtain the next iteration of the Markov 
chain, a sample is drawn from a proposal distribution, and it is accepted or 
rejected according to a certain probability so that the target distribution is 
preserved. In Gibbs sampling, the parameter 9 is partitioned into several com- 
ponents, and at each iteration, we update each component in turn by drawing 
from its conditional posterior distribution given all other components. The actual 
MCMC algorithm used for fitting our two-stage model is a hybrid algorithm 
that cycles through the coordinates of the parameter vector in a Gibbs sampling 
fashion but uses an M-H strategy for each conditional draw. The algorithm is 
carefully monitored; several Markov chains from different starting values are run 
to ensure that they reach the same target distribution. 
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3. Joint Fitting: Results and Discussion 

We now discuss the results from the fitted hierarchical model. This is a more 
coherent analysis than the separate fitting of the two stages as described in 
Section |2~T1 Thus, the results in this section represent our final estimates. 

3.f . Cyclc-to-Cyclc Dependencies 

Of particular interest are the estimates of the three second-stage relationships. 
These estimates are summarized in Table H3 based on an MCMC calculation 
that fits the two stages jointly, using available data from January 1749 to Oc- 
tober 2011. As noted before, because of the large number of observations, cycle 
characteristics such as Ci,tg ,£mLc, and if are well constrained by the first- 
stage model (Equation ) , and hence the estimates reported in Table [5] are 
quite close to those based on an ordinary least squares (OLS) regression using 
fixed Ci,tg , £max>ii • For example, with a, , (tL fixed at their respective 
posterior means, the OLS estimates (standard errors) of 62 and 72 are 8.5(±0.8) 
and — 0.43(±0.08), respectively. The standard errors in Tableware slightly larger 
because they account for the extra uncertainty in estimating Cj, , imax, and . 
The quality of the predictive relationships is shown in Figure [SJ where we display 
the ratios of the values predicted from the previous cycle to those measured for 
that cycle, together with 1-cr prediction error bars. Here the posterior mean 
estimates (as in Table [5]) are regarded as measured or estimated values; the 
predicted ones are computed using the equations in Tabled The prediction error 
bars are computed by simulation. Note that the coefficients Si, 7$, i = 1,2,3, 
are estimated using all cycles (including those that come after the one being 
predicted). 

Table 6. Summary of three relationships and their parameter estimates when fitting 
the two stages jointly. The fitted values are posterior means and the standard errors 
are posterior standard deviations. 



Cycle parameter 


Relationship 




Fitted value (Std. err.) 


Amplitude 


Cj+l ~ 5i + 7i c i/(*o 


hi) s 
L max } 


<5i = 4.1 (±1.5) 
71 = 3.9 (±1.0) 


Time to maximum 


tmax — *Q ~ fl 2 + 


■72C i+ i 


<5 2 = 8.5 (±1.0) 

72 = -0.43 (±0.09) 


Time to minimum 


t (i+l) ^,5,1 


■73C i+ i 


<5 3 =4.3 (±1.5) 
73 = 0.31 (±0.15) 



3.2. Predictions for Cycle 24 

One advantage of a Bayesian hierarchical model in this context is that, once 
we obtain the samples from the posterior distribution, prediction of the char- 
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acteristics of the current (incomplete) cycle is obtained automatically. Table [7] 
reports summaries of the posterior inference for Cycle 24, using data up to 
November 2008, May 2010, and October 2011, respectively. Based on data up to 
October 2011, Cycle 24 is estimated to rise to maximum in January - February 
2014 ± six months, with a maximum smoothed monthly sunspot number of 
97 ± 15, where the estimates are posterior means and the error bars are posterior 
standard deviations. (The maximum smoothed SSN, or the expected SSN at solar 
maximum, is (f3 + Ci) 2 + <7 2 , after accounting for the square-root transformation 
in Equation ([3]).) It is likely to be a weak cycle with a longer-than- usual (an 
expected 12.1 years) total duration, although the uncertainty associated with 
this estimate is fairly large. We observe that the estimated maximum smoothed 
sunspot number is relatively stable across the three analyses. The large error bars 
associated with the November 2008 analysis highlight the inherent difficulty in 
making predictions before or at the onset of a cycle. 

Table 7. Cycle 24 predictions based on MCMC fitting of the two stages jointly. Fitted 
values are posterior means and standard errors are posterior standard deviations. Max. SSN 
refers to the smoothed monthly average sunspot number at the peak of Cycle 24. Time to 

rise [years] is defined as tmax — *q 24 ' ■ 





C24 


Max. SSN 


Time of max. [yrs] 


Time to rise 


Cycle length 


Nov 08 


9.0 ±1.6 


96 ±32 


Mar 2013 ± 0.98 


4.7 ± 1.1 


11.9 ± 1.7 


May 10 


8.2 ± 1.2 


77 ± 21 


May 2014 ± 0.69 


5.3 ±0.84 


12.1 ± 1.6 


Oct 11 


9.3 ±0.77 


97 ± 15 


Jan/Feb 2014 ± 0.48 


4.8 ±0.55 


12.1 ± 1.5 



Figure [6] illustrates the estimates of the averages of Y t in Equation <j3j> (specif- 
ically, (/3 + Ut) 2 + <J 2 )- The solid curve represents the posterior mean while the 
upper (lower) dashed curve represents the 95% (5%) posterior quantile. Note 
that estimates for time points in the past are well constrained because of the 
available data, but future predictions are much more variable. The two-stage 
model is well suited for combining two pieces of information that have potential 
predictive power: sunspot number observations that clearly belong to Cycle 24, 
and the prescription of the second-stage model (Equations (HJ) - ©) which relates 
the characteristics of Cycle 24 to those of previous cycles. Given the relatively few 
observations at the beginning of Cycle 24, the predictions arc heavily influenced 
by the second-stage model. As Cycle 24 progresses, direct observations will play 
a heavier role and the uncertainties associated with the predictions will diminish. 
This is illustrated by the reduction in the uncertainty band in the bottom panel 
which includes 35 more months of observations (up to October 2011) compared 
to that in the top panel (up to November 2008). (This reduction in uncertainty 
is also apparent from Table UJ) When the more recent data are included, the 
predictions are more driven by direct observations from Cycle 24; the fitted 
values are similar, but the 90% predictive intervals are appreciably narrower. 
This shows that the latest data are reasonably consistent with the second-stage 
relationships, and combining the two stages shrinks the error bars. 
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4. Summary 

Wc have carried out a comprehensive statistical analysis of the sunspot record. 
After suitably transforming the data to stabilize variance, we parameterized the 
shape of a cycle by its amplitude (maximum in the sunspot number), time to 
rise to maximum, time to fall to minimum, and the gap between its end and 
the start of the next cycle. By computing correlations between these parameters 
both within each cycle and between adjacent cycles, we have derived a set of 
three predictive relations. These relations are ordered, i.e. sequential: amplitude 
must be predicted first before duration and rise time. Correlations that depend 
on computing amplitude second are not robust, as they are subject to two 
influential points from early in the sunspot record (see also Vaquero and Trigo 
2008). Analyses carried out in a different order will thus lead to spurious results. 

These relations can be used to predict the values of the parameters for the 
following cycle. We find that the best estimate for the peak in Cycle 24 is in 
early 2014, with an uncertainty of half a year. The maximum in the smoothed 
sunspot number record is expected to be ps 97 ± 15, and the cycle is expected 
to last 12.1±1.5 years from the latest solar minimum (approximately November 
2008). These are in the middle of the range of predictions in the literature prior 
to the onset of Cycle 24, and are consistent with the current estimates of the 
cycle parameters. 

Wc have searched for, but do not find, any evidence for persistence beyond one 
cycle. There is no predictive power beyond the cycle that follows; no correlations 
are present, and the cycles do not retain any memory. 

We also find that the cycles do not ever vanish completely. We find statistical 
evidence that the next cycle usually begins before the current cycle ends, as the 
gap between the end of a cycle and the start of a new one is usually negative. 
Furthermore, we find that sunspots do not vanish entirely even if the gap were 
positive; however, the data are not sufficient to tell whether this holds true even 
in the absence of activity cycles. 
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Figure 5. Ratios of predicted parameter values to measured values. The values predicted for 
the parameters in each cycle (amplitude, rise time and fall time) based on the parameters of 
the previous cycle, are compared with the values directly estimated for that cycle. 
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Figure 6. Predictions of Cycle 24 obtained by fitting the two stages jointly, using data up to 
November 2008 (top), May 2010 (middle) and October 2011 (bottom). The posterior mean of 
monthly average sunspot numbers is shown as the solid curve, and the 5% and 95% posterior 
quantiles are shown as dashed curves. The top figure illustrates predictions for a completely 
new cycle, and the bottom for one that is well in progress. May 2010 is chosen as it lies half 
way in between November 2008 and October 2011. Note that the uncertainty in the predictions 
is reduced considerably when more data from the current cycle are included. 
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